SUMMARY 


In  this  final  report  covering  the  period  1993-1994,  the  preliminary  results  of  the 
research  on  direct  numerical  simulation  of  compressible  turbulent  flows  are  presented. 
Specifically,  the  adaptive  Legendre  polynomial  spectral  element  method  with  the  mked 
explicit-implicit(MEI)  Taylor-Galerkin  formulation  is  used  in  order  to  capture  detailed 
physics  involved  in  three-dimensional  shock  wave  turbulent  boundary  layer  interactions. 

The  main  activity  for  the  second  year  was  concerned  with  construction  of  three- 
dimensional  data  structure  of  the  computer  program.  The  material  covered  in  this  final 
report  consists  of  (1)  variations  of  fluctuation  velocity  components,  Reynolds  stress 
contours,  and  turbulent  kinetic  energy  variations,  corresponding  to  the  two-dimensional 
shock  wave  turbulent  boundary  layer  interactions  examined  in  the  first  year  progress 
report  and  (2)  the  preliminary  results  of  the  three-dimensional  flow  including  the  vorticity 
distributions.  Studies  on  three-dimensional  turbulent  statistics  such  as  power  spectral 
densities  vs  wave  numbers  (or  fi'equencies),  spectral  energy  transfer  and  cascade,  the 
effects  of  production  and  dissipation,  inertial  subrange,  etc.  are  currently  in  progress  and 

will  be  included  in  the  next  progress  report. 

In  conclusion,  the  proposed  direct  numerical  simulation  of  compressible  turbulent 
flows  using  the  adaptive  Legendre  polynomial  spectral  element  with  the  MEI  Taylor- 
Galerkin  formulation  has  been  shown  to  be  robust,  accurate,  and  efficient.  Further 
research  is  required,  however,  to  demonstrate  its  ultimate  impact  on  CFD  through 
exhaustive  examples  and  comparisons  with  available  experimental  data  in  the  future. 
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I .  INTRODUCTION 

The  last  decade  has  seen  unprecedented  technological 
innovations  in  computational  fluid  dynamics,  prompted  partic¬ 
ularly  by  the  increase  in  technical  requirements  of  aerospace 
research.  Namely,  the  flowfields  due  to  high  velocities, 
compressibility,  shock  waves,  turbulence,  and  high  temperature 
have  been  the  focus  of  intensive  research  in  the  past  [1-9]. 

When  shock  waves  interact  with  turbulent  boundary  layers  in 
external  or  internal  flows  special  considerations  are  required 
due  to  widely  disparate  time  and  length  scales,  corresponding  to 
different  physical  phenomena  -  namely,  turbulence  microscales  and 
shock  wave  surface  discontinuities.  Here  we  are  faced  with  the 
smallest  time  and  length  scales  which  may  severely  affect  the 
computational  requirements.  To  cope  with  such  requirements 
various  numerical  strategies  have  been  developed  using  finite 
difference  methods  (FDM)  and  finite  element  methods  (FEM) . 
Incorporated  or  implemented  into  either  FDM  or  FEM  are  the  finite 
volume  methods  (FVM)  and  spectral  element  methods  (SEM) . 

Modeling  of  turbulence  has  been  the  controversial  subject. 

Closure  models,  probability  density  functions  (PDF),  large  eddy 
simulation  (LES) ,  direct  numerical  simulation  (DNS) ,  and  other 
methods  have  been  reported. 

The  purpose  of  the  present  study  is  to  demonstrate  the 
superiority  of  DNS  combined  with  unstructured  adaptive  spectral 
element  methods  in  dealing  with  combined  turbulence  and  shock 
waves  for  both  internal  and  external  flows  of  aerospace  vehicles. 
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This  work  is  motivated  by  the  fact  that  DNS  can  be  achieved  via 
adaptive  h-p  methods,  combining  the  mesh  refinement  (h-method) 
with  spectral  polynomial  degrees  of  freedom  (p-method) .  It  is 
well  known  that  the  most  crucial  aspect  of  turbulent  flows  is 
microscales  involved  in  boundary  layers  (viscous  sublayer,  buffer 
zone,  and  turbulent  core) .  This  is  where  the  spectral  polynomial 
degrees  of  freedom  can  be  increased  as  desired  since  the  mesh 
refinement  alone  is  incapable  of  resolving  the  microscale 
requirements.  In  this  way,  turbulence  modeling  techniques  can  be 
avoided.  Furthermore,  the  current  practice  in  DNS  to  use 
extensive  refinements  in  finite  difference  discretization  may 
also  be  avoided.  Babuska  and  his  co-workers  [10-12]  and  Oden  and 
his  co-workers  [13-15]  contributed  to  the  advancement  of  FEM  h-p 
adaptive  methods.  Their  applications  have  not  been  extended  to 
shock  waves  interacting  with  turbulent  boundary  layers. 

Chung,  et  al.  [16-19]  have  studied  finite  element  strategies 
as  applied  to  shock  wave  turbulent  boundary  layer  interactions  in 
reacting  flows  and  explored  applications  of  direct  numerical 
simulation  in  characterizing  the  shock  wave  turbulent  boundary 
layer  interaction.  The  main  emphasis  in  the  present  study  is  to 
establish  the  basic  theory  and  computational  strategies  involved 
in  the  Legendre  polynomial  spectral  element  method  and  to  present 
preliminary  computational  results.  Development  of  theory  and 
formulations  include  irregular  node  connectivity  of  Legendre 
polynomials  of  various  orders.  Comparisons  with  experimental  re¬ 
sults  have  demonstrated  superiority  of  the  direct  numerical  simu- 
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lation  over  the  standard  K-e  model  with  compressibility  effects. 

In  what  follows  we  discuss  the  governing  equations  and 
solutions  of  Navier-Stokes  system  of  equations  via  mixed 
explicit-implicit  (MEI)  Taylor-Galerkin  methods  in  Section  2, 
direct  numerical  simulation  with  h-p  adaptive  methods  in  Section 
3 ,  calculations  of  DNS  perturbation  variables  in  Section  4 , 
calculations  of  f lowf ield-dependent  implicitness  parameters  and 
stability  analysis  in  Section  5,  numerical  applications  in 
Section  6,  and  conclusions  in  Section  7. 


II.  GOVERNING  EQUATIONS  AND  SOLUTIONS  OP  NAVIER-STOKES 
SYSTEM  OF  EQUATIONS 

A  convenient  form  of  governing  equations  for  compressible 
viscous  flows  may  be  written  in  terms  of  conservation  variables 
as  follows: 

(1) 
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The  solution  of  governing  equations  will  be  carried  out  using  the 
Taylor-Galerkin  approach  with  test  and  trial  functions  given  by 
isoparametric  and  Legendre  polynomials  by  means  of  mixed 
explicit/ implicit  schemes.  In  general,  explicit  schemes  are 
inexpensive  but  less  accurate  in  comparison  with  implicit  schemes 
for  regions  of  high  pressure  or  velocity  gradients.  In  case  of 
rapid  variations  of  gradients  throughout  the  domain  it  is  often 
desirable  to  devise  a  scheme  in  which  implicitness  can  be 
adjustable  in  accordance  with  gradients,  more  implicit  for  the 
region  of  high  gradients  and  less  implicit  or  fully  explicit  for 
the  region  of  low  gradients. 

In  expanding  in  Taylor  series  about  U",  we  introduce  the 
implicitness  parameters  s^  and  Sj  for  the  first  and  second 
derivatives  of  U  with  respect  to  time  [17-19],  respectively, 


+  At 


dt 
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Substituting  (3)  into  (2)  yields 
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It  follows  from  (1)  that 
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dU 

dt 


dr, 


dG, 


dxi  dxj 


+  B 
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Here  P.  is  a  function  of  U  and  G,-  is  a  function  of  U  and  its 
gradient  U . ,  so  that  we  denote  the  convective  Jacobian  a-, 
dissipative  Jacobian  bj,  and  dissipative  gradient  Jacobian  as 
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The  second  derivative  of  U  with  respect  to  time  may  now  be 
written  in  terms  of  these  Jacobians. 
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Substituting  (5)  and  (6)  into  (4)  yields 
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In  order  to  provide  different  implicitness  (different  amount 
of  damping  or  dissipation)  to  different  physical  quantities,  we 
reassign  s^  and  Sj  associated  with  G-,  respectively 

s^AG^  -  ,  SjAG^  -  s^AG^ 

with  the  various  implicitness  parameters  defined  as 

Si  =  first  order  convective  implicitness  parameter 
Sj  =  second  order  convective  implicitness  parameter 
S3  ~  first  order  dissipative  implicitness  parameter 
s^  =  second  order  dissipative  implicitness  parameter 

These  implicitness  parameters  play  a  significant  role  with 


(8a, b) 


(8c) 
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s,  and  Sj  controlling  shock  discontinuity  resolutions  and  S3  and 
s^  dictating  turbulent  eddy  microscale  resolutions.  The 
implicitness  parameters  are  considered  to  be  dependent  not  only 
on  time  but  also  on  space.  In  particular,  s^  and  Sj  are 
associated  with  Mach  number  changes,  whereas  S3  and  s^  depend  on 
Reynolds  number  changes  between  time  steps  and  between  upstream 
and  downstream  locations.  Further  details  of  these  implicitness 
parameters  are  given  in  Section  5. 

Neglecting  the  third  order  spatial  derivatives  of 
conservation  variables  associated  with  c,-|j  in  (7)  and 
substituting  (8a,  b)  into  (7)  lead  to  the  residual. 
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(9) 


where  all  Jacobians  a,.,  bj  and  Cjj  are  assumed  to  remain  constant 
spatially  within  each  time  step  and  to  be  updated  at  subsequent 
time  steps. 

The  Galerkin  integral  of  (9)  may  now  be  carried  out  as 
follows: 

f$„lZ(Z7,  F^,  a^)dCl  =  0 
0 


(10) 
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where  refers  to  the  spectral  test  functions,  and  the 
conservation  variables  are  interpolated  by  the  spectral  trial 
functions  given  by  the  combination  of  isoparametric  functions 
and  Legendre  polynomials  as 

t7(x,  t)  =  *.  (x)  IT,  ( t) ,  Fiix.t)  =  .  Oi(x,  t)  =  *.(x)o.i(t)  (11) 


Substituting  (11)  into  (9)  and  (10)  yields 
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Here  terms  on  RHS  of  (12)  are  assembled  into  boundary  nodes 
of  B  .  on  LHS.  Notice  also  that  indices  i,  j,  k  =  1,  2 

AdPS 
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associated  with  the  Jacobians  imply  directional  identification  of 
each  Jacobian  matrix  (a,,  aj,  b, ,  h^,  c,2/  ^21'  ^22^  with  r,s  = 

1,  2,  3,  4  denoting  entries  of  each  of  the  4x4  Jacobian  matrices. 
It  should  be  warned,  however,  that  these  Jacobian  matrices  must 
be  multiplied  precisely  as  dictated  by  summing  through  repeated 
indices,  not  through  matrix  multiplications  as  a  whole. 

It  is  interesting  to  note  that  all  implicitness  parameters 
can  be  shown  to  be  functions  of  flowfields  between  upstream  and 
downstream,  and  that  the  convection  implicitness  parameters  s^ 
and  Sj  associated  with  the  first  term  in  are  analogous  to 

the  total  variation  diminishing  (TVD)  limiter  in  the  FDM 
literature  (see  Appendix  A) .  With  an  adequate  choice  of  these 
implicitness  parameters,  acceptable  resolutions  of  shock  waves 
have  been  verified. 

On  the  other  hand,  the  diffusion  implicitness  parameters  S3 
and  s^  are  capable  of  alleviating  and  accommodating  the  stiffness 
involved  in  turbulent  diffusion  or  finite  rate  chemistry  (if 
reacting) .  No  analogy  can  be  shown  since  they  do  not  exist  in 
other  numerical  schemes.  It  should  also  be  noted  that 
interactions  between  convection  and  diffusion  are  achieved  by 
means  of  the  terms  associated  with  the  products  aj^^^  bj^^^  and  b.^^ 
ajgq.  These  terms  are  particularly  important  for  shock  wave 
turbulent  boundary  layer  interactions  where  the  effect  of 
convection  upon  diffusion  and  vice  versa  is  crucial  in  order  to 
resolve  turbulence  microscales  as  disturbed  by  shock  wave 
interactions.  We  shall  refer  to  these  terms  as  convection- 
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diffusion  interaction  terms. 

For  low  Mach  nximber  flows  density  variations  become 
insignificant.  In  this  case  the  flow  is  essentially 
incompressible  and  checkerboard  type  pressure  oscillations  occur 
unless  a  suitable  computational  scheme  is  devised  to  ensure  the 
mass  conservation.  In  the  present  study  the  implicitness 
parameters  as  calculated  from  the  flowfield  are  expected  to 
accommodate  the  conservation  of  mass  and  prevent  pressure 
oscillations.  In  this  vein,  the  effect  of  compressibility  will 
also  be  automatically  accommodated  for  high  Mach  number  flows. 

III.  DNS  VIA  UNSTRUCTURED  h-p  ADAPTIVE  METHODS 

Our  objective  here  is  to  resolve  time  and  length  scales 
involved  in  turbulence  interacting  with  shock  waves  using  the 
unstructured  h-p  adaptive  spectral  element  methods.  One  approach 
is  to  refine  the  mesh  (h-methods)  until  further  refinement  is 
unproductive,  at  which  time  the  spectral  degrees  of  freedom  (p- 
methods)  are  increased  in  order  to  reduce  errors  as  desired,  such 
as  in  the  region  of  turbulent  viscous  sublayer.  However,  the 
more  desirable  approach  is  to  optimize  between  the  mesh 
refinement  and  spectral  orders.  Thus,  the  most  crucial  aspect  of 
the  h-p  methods  is  to  determine  the  best  possible  change  in  the 
mesh  structure  to  reduce  the  local  error  to  a  minimum.  Should  h 
(mesh  size)  be  decreased  or  should  p  (polynomial  or  spectral 
degrees  of  freedom)  be  increased?  Although  some  work  in 
optimization  between  the  h-  and  p-  processes  has  been  reported 
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[13-15],  the  subject  of  optimization  appears  to  be  an  open 
question.  Thus,  our  approach  in  this  study  is  to  refine  the  mesh 
until  shock  waves  are  adequately  computed  and  then  resort  to  the 
p-  version  with  higher  order  Legendre  polynomials  in  order  to 
resolve  turbulence  microscales.  Toward  this  end,  the  error 
indicator  e  may  be  defined  in  terms  of  density  for  shock  waves 
and  velocity  for  turbulence.  Some  of  the  options  are  given  as 
follows: 

0  =  h  "o^lpl^a/lpl^i 


0  =  h  Q  \Vi\g2/\vi\if. 


(13b) 


where  h  is  the  mesh  parameter  and  various  Sobolev  space  (H^) 
seminorms  are  defined  as 


(14b) 


The  choice  among  these  options  depends  on  various  physical 


aspects  of  the  given  problem,  whether  local  errors  are  dominated 
by  density,  velocity  components,  their  gradients,  or  second 
derivatives.  For  the  purpose  of  the  examples  dealt  with  in 
Section  5,  we  utilize  Eg.  (13a)  and  Eq.  (13b)  for  the  h- 
adaptivity  associated  with  shock  waves  and  the  p-  adaptivity 
associated  with  turbulent  boundary  layers,  respectively. 

Direct  numerical  simulations  for  turbulent  flows  are 
achieved  by  higher  spectral  orders  using  Legendre  polynomials 
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[10-12].  The  spectral  element  interpolation  functions  for  the 
corner  nodes,  side  modes,  .and  interior  modes  in  a  two-dimensional 
4-node  isoparametric  element  are  as  follows  (see  Fig.  la) : 

Four  corner  node  functions: 

(i-Ti)  , 

^  "  (15) 

(1+Tl)  , 

4  4 


Side  mode  functions: 
jg(SM)  =  4(P-1)  ,  i=2,  . 


*i"’  =  |(l-tl)G„(5) 


•  p 

*(S2) 


A(l-5)G,„(tl) 


(16) 


Interior  mode  functions  (bubble  functions) : 

j,(l)  .  (P-2)  (gjl.)  _ 

2 

'  m.n  =  2....p,  m+n^p 


(17) 


where  Ni?’’  and  denote,  respectively,  the  total  number  of 
interpolation  functions  available  for  sides  n  (1,2, 3, 4)  and  the 
interior,  p  is  the  highest  Legendre  polynomial  order  chosen,  and 
the  Legendre  interpolation  functions,  G^($)  and  G^(ti),  are 
defined  in  terms  of  the  Legendre  polynomials  p,(5)  as  follows; 


1 

y/2  {201-1) 


(Pn.(5)-P.-2(0) 


(18) 


The  same  procedure  applies  to  the  n -direction. 

Notice  that  the  corner  node  interpolation  functions  are 
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linear  whereas  the  side  inodes  and  interior  inodes  depend  on  the 
order  p  of  Legendre  polynomials  chosen.  For  a  three-dimensional 
domain  tensor  products  in  the  third  direction  C  for  a  hexahedral 
element  can  be  constructed  similarly  as  in  the  two-dimensional 
case,  but  with  surface  mode  functions  as  well  as  edge  and 
interior  mode  functions. 

As  a  conseguence  of  these  Legendre  polynomials  representing 
the  side  and  interior  mode  functions,  any  variable  U  may  be 
modeled  as 

u  =  (19) 

where  *  and  O”  denote  the  side  and  interior  mode  functions, 

Ri  inn 

respectively,  and  0^  and  0^  are  spectral  coefficients  to  be 
calculated  from  additional  integrals, 

=  0  (20a) 

f  <5^’jZ(U)dQ  =  0  (20b) 

JQ 

Here  the  standard  static  condensation  may  be  used  to  eliminate  0^ 
and  U  .  Thus  the  final  form  of  the  finite  element  equations  is 

riffi 

modified  to 


(A.p6rp«  +  ^rprs)  ^  '  <1:  +  +  ^ar  ^ 

where  W  "  acts  as  a  source  term  reflecting  the  Legendre 

<zr 

polynomial  side  and  interior  mode  functions  obtained  from  the 
static  condensation  of  (12),  (20a),  and  (20b).  In  this  way, 
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variables  are  calculated  only  at  the  corner  nodes ,  regardless  of 
high  order  Legendre  polynomial  approximations. 

Our  objective  here  is  to  satisfactorily  simulate  turbulent 
microscales  within  an  element.  All  side  mode  and  interior  mode 
interpolation  functions  vanish  at  the  corner  nodes  but  exhibit 
high  frequency  variations  according  to  the  order  of  Legendre 
polynomials  along  the  sides  and  interior  domain.  It  is  intended 
that  such  Legendre  polynomial  microscales  be  capable  of 
simulating  the  physical  microscales  of  turbulence  such  as  those 
of  Kolmogorov  and  Taylor  (Fig.  lb) .  Microscales  involved  in 
viscous  sublayer,  buffer  zone,  and  turbulent  core  are  of  the 
order  of  lO"’  mm  whereas  characteristic  lengths  of  domain  of 
interest  may  be  over  meters  or  kilometers.  Thus  the  h-adaptivity 
alone  is  severely  limited  and  naturally  we  seek  a  remedy  of  this 
situation  in  the  h-p  adaptivity  utilizing  the  highest  spectral 
orders  required  for  accuracy. 

Transition  elements  involve  irregular  nodes  in  the  h- 
refinement  process  —  node  c  for  element  T-1  and  nodes  c  and  d 
for  element  T-2  as  shown  in  Fig.  Ic.  To  assure  linear 
approximations  for  both  unrefined  and  refined  elements,  we  must 
eliminate  irregular  nodes  involved  in  the  unrefined  elements. 

This  can  be  done  simply  by  assuming  [13] 

(22) 

This  will  lead  to 


{N,M  =  1,2, 3, 4) 


(23a) 
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U^r-2)  ^  h^-2)  U^T-2)  {n.M  =  1,2,3,A)  (23b) 

where  ’’  and  denote  the  auxiliary  matrices  for  T-1  and 

T-2  elements,  respectively^ 
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[2  2  J 

and  and  are  the  unknown  nodal  vectors  in  terms  of  the 

global  nodes  Q,  R,  and  S  with  the  irregular  nodes  c  and  d 
eliminated. 

The  advantage  of  Legendre  polynomials  is  an  ease  in  dealing 
with  side  and  interior  modes  which  do  not  require  specification 
of  nodes  physically  located.  This  is  especially  beneficial  for 
side  modes  in  establishing  boundary  continuities.  Continuity  of 
variables  and  gradients  along  the  inter-element  boundaries  is  to 
be  dictated  by  the  higher  order  polynomials  between  the  two 
adjacent  elements. 


Similar  extensions  can  be  made  for  three-dimentional  geome¬ 
tries.  The  three-dimensional  Legendre  polynomial  spectral  element 
configuration  is. .shown  in  Fig.  Id.  It  is  seen  that  side  modes 
and  the  interior  mode  of  a  2-D  element  become  edge  modes  and 
face  modes  for  a  3-D  element,  respectively.  These  modes  are  in 
addition  to  the  3-D  interior  mode. 


IV.  DNS  PERTURBATION  VARIABLES 
It  is  well  known  that  DNS  is  expected  to  provide  information 
in  turbulence  microscale  levels  at  the  expense  of  excessive 
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refinements  of  domain  discretization  [6].  The  purpose  of  the 

present  study  is,  instead,  to  avoid  such  refinements  by  means  of 

implementing  high  spectral  Legendre  polynomial  orders.  The 

Navier-Stokes  solver  as  introduced  here  allows  unsteady  time 

accurate  solutions  from  which  perturbation  variables  (f')  can  be 

calculated  as  the  difference  between  the  Navier-Stokes  solution 

(f)  and  its  time  average  f  [20,  21], 

f  =  f-f  (25) 

This  computation  can  be  conducted  throughout  the  Navier-Stokes 
integration  time  steps  or  upon  arrival  at  quasi-steady  state. 
Strictly  speaking,  in  shock  wave /turbulent  boundary  layer 
interactions  a  complete  steady  state  is  never  realized  as 
unsteady  eddy  motions  persist  indefinitely,  although  background 
flowfields  may  become  steady.  This  is  referred  to  as  the  quasi¬ 
steady  state.  The  time  average  of  Navier-Stokes  solution  is 
performed  using  the  Gaussian  quadrature.  In  this  process 
complicated  physical  phenomena  such  as  homogeneous  and 
inhomogeneous,  isotropic  and  anisotropic,  and  nonstationary 
nature  of  perturbation  flowfields  for  a  compression  corner  with 
shock  wave  turbulent  boundary  layer  interactions  can  be  resolved. 

Furthermore,  all  perturbation  variables  as  calculated  from 
(25)  can  be  transformed  via  fast  Fourier  transform  to  generate 
power  spectral  density  vs  frequency  domain.  In  section  6  we 
examine  various  perturbation  variables  as  well  as  background 
flowfield  data.  As  a  result  of  this  study,  more  details  of  shock 
wave  turbulent  boundary  layer  interactions  such  as  variations  of 
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turbulent  kinetic  energy  vs  shock  strength,  laminar-turbulence 
transition  instability,  relaminarization,  effects  of  dilatation, 
etc. ,  can  be  rigorously  examined  in  comparison  with  the  previous 
investigations  [6,  22-23].  Some  limited  results  and  discussion 
of  these  subjects  are  presented  in  Section  5. 


V.  CALCULATIONS  OF  FLOWPIELD-DEPENDENT  IMPLICITNESS 
PARAMETERS  AMD  STABILITY  ANALYSIS 

The  success  of  hp  version  spectral  element  method  described 
above  depends  on  accurate  calculations  of  flowfield— dependent 
implicitness  parameters.  With  appropriate  choices  of  nodal 
points  in  a  one-dimensional  case  it  was  demonstrated  in  Section  2 
that  the  terms  associated  with  the  convection-implicitness 
parameters  (s,,  Sg)  in  the  MEI  formulation  are  analogous  to  the 
the  FDM-TVD  methods.  To  further  examine  both  convection  and 
diffusion  implicitness  parameters  we  propose  the  following 
criteria: 


Si,  Sj 


min(r, 1)  r  >  a 
0  r  <  a 

1  ^  =  0 


(25) 


with 


r  = 


AAf 


(26) 


where  AM  is  the  difference  between  the  maximum  and  minimum  Mach 
number  (AM  =  M  -  M  ,  )  within  a  finite  element,  and  a  is  a  user- 


specified  small  number  (a  -  0.01). 
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f  min(s,  1)  s  >  P 
s,.s,  =  0  s  <  p 

[  1  ■^®(ndn)  “  ® 


(27) 


with 

s  =  p=0.01  (28) 

■^®(inin) 

where  ARe  is  the  difference  between  the  maximum  and  minimum 
Reynolds  number  (ARe^^,,  -  ARe^^,„,)  within  a  finite  element  and  6 
is  a  user-specified  small  number  (6  -  0.01).  For  reacting  flows, 
the  diffusion  parameters  S3  and  s^  are  also  calculated  from 
Damkohler  nximbers  using  the  similar  criteria  as  in  (27)  and  (28)  . 
Depending  upon  the  length  and  time  scales  involved  in  the  actual 
flowfield,  the  diffusion  implicitness  parameters  are  governed  by 
either  turbulence  or  finite  rate  chemistry,  whichever  are  larger 

if  chemical  reactions  are  included. 

As  a  result  of  both  convection  and  diffusion  implicitness 
parameters  based  on  the  flowfield  variables,  it  is  possible  to 
derive  an  expression  for  the  amplification  factor  for  (12)  in  one 
dimension  in  the  form 

l^|(l-s,)C^^2(l-s.)-^-3(l-5«)  ^(C03»-1)  ^  J(s,-l)C3in4. 

+  JSiCsin* 

[  ^  Re  Re^ 

with 

Be  =-5*, 

h  V 

For  illustration,  the  amplification  factors  for  various 
courant  number  C  are  shown  in  Fig.  2.  It  is 


modes  m  versus 
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evident  that  stability  conditions  are  much  more  favorable  than  those  given  in  [32]  where 
only  the  convection  implicitness  parameters  were  considered.  It  is  also  shown  that  as  the 
Reynolds  number  decreases  stability  increases  for  intermediate  modes  (m=0.5)  with  an 
increase  of  difiusion  implicitness  parameters  s3  and  s4. 

VL  APPLICATIONS 

In  the  1993  progress  report,  the  validity  of  the  use  of  Legendre  polynomial 
spectral  elements  in  conjunction  with  adaptive  mesh  refinements  was  demonstrated  using 
two  example  problems  (flat  plate  and  compression  comer).  Fourth  order  Legendre 
polynomials  were  used  for  direct  numerical  simulation  of  these  problems.  In  that  report 
the  results  of  turbulent  statistics  were  not  presented. 

In  Fig.  3,  the  perturbution  velocity  components  corresponding  to  the  flate  plate 
shock  wave  turbulent  boundary  interactions,  page  37,  the  1993  progress  report  are  shown. 
These  variables  are  in  a  random  unsteady  state  within  the  boundary  layer  (solid  line)  and 
more  or  less  steady  state  outside  the  boundary  layer  (dashed  line).  Here,  one  unit  of 
nondimensional  time  is  equivalent  to  0.96  ms,  and  thus  the  unsteady  motion  frequencies 
vary  from  approximately  62  to  135  kHz  for  the  streamwise  component  and  31  to  94  kHz 
for  the  spanwise  component. 

Fig.4  shows  perturbation  normal  and  shear  stress  components,  indicating  that 
shock  wave  turbulent  boundary  layer  interactions  contribute  to  dilation  (compressibility 
effect)  as  evidenced  by  greater  normal  stresses  than  sheal  stresses. 

Also  of  interest,  as  shown  in  Fig.  5,  is  the  sharp  rise  of  turbulent  kinetic  energy  in 
the  boundary  layer,  away  from  the  wall  but  slightly  below  the  boundary  layer,  where 
shock-turbulence  interactions  are  maximum,  indicating  the  amplification  of  turbulent 
kinetic  energy. 
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The  major  achievement  of  the  past  year  is  the  completion  of  computer  algorithm 
for  the  three-dimensional  data  structure.  To  check  out  this  algorithm  a  three-dimensional 
sharp-leading-edged  fin  has  been  investigated  for  swept  shock  wave  turbulent  boundary 
layer  interactions.  Fig.6a  shows  the  physical  domain  for  a  3-D  sharp  fin  (  a=20  )  and  a 
general  flowfield  structure  (Fig.6b,c)  [33].  The  inlet  boundary  conditions  are  the  same  as 
used  by  Knight  et  al.  [34]  and  the  corresponding  flowfield  structure.  Here  the  free  stream 
Mach  number  and  temperature  are  M*=2.93  and  T*  =92.39  K,  corresponding  to  the 
chamber  pressure  and  temperature  of  680  kPa  and  25 1  K,  respectively,  with  the  Reynolds 
number  being  7  x  lo*  /m.  The  boundary  layer  thickness  Sg  at  the  apex  of  the  fin  is  1.4cm, 
yiftlHing  a  Reynolds  number  Re^g=  9.8  x  10^.  In  order  to  match  the  boundary  conditions  as 

used  for  the  experiments  [34]  the  flowfied  behind  the  fin  is  calculated  as  a  flat  plate 
boundary  layer  such  that  the  computed  boundary  layer  thickness  Sg  is  set  equal  to  the 
experimental  value  of  1.4  cm.  On  the  solid  surfaces  no  slip  and  adiabatic  wall  boundary 
conditions  are  applied.  On  the  upper,  lateral,  and  at  downstream  exit  boundaries  the  flow 
variables  are  set  free.  Adaptively  spaced  grid  points  are  33,  41,  and  31  in  the  streamwise, 
spanwise,  and  vertical  directions,  respectively.  No  attempt  is  made  for  further  adaptive 
refinements  for  the  geometry  shown  in  Fig.6d  at  this  time. 

Fig.  7  shows  the  background  flowfield  based  on  the  geometric  configurations  and 
boundary  conditions  described  in  Fig.6,  as  observed  from  the  front  (x-z  and  y-z  faces).  As 
such,  no  details  of  the  hidden  portion  are  shown.  Subject  to  the  detailed  examination  in 
comparison  with  the  experimental  data  to  be  given  in  Fig.  8,  we  are  only  able  to  detect  a 
general  qualitative  trend  of  the  flowfield  in  Fig.7.  It  is  to  be  noticed  that  the  trend  is  in 
reasonable  agreement  with  the  results  of  Narayanswami  et  al.  [35],  with  density  and 
pressure  increasing  drastically  along  the  shock  waves,  the  temperature  rise  being 
distributed  along  the  flat  plate,  and  Mach  number  sharply  decreasing  through  the  shock 
waves  toward  the  flat-plate  boundaiy. 
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Fig.  8  shows  comparisons  of  yaw  angles  at  various  locations  :  (a)  x=7.549  cm, 
y=-6. 1087cm  and  (b)  x=8.827cm,  y=-7.3787cm  with  the  experimental  and  other 
computational  data  [34].  As  pointed  out  in  [34]  uncertainties  in  the  experimental  data 
(portion  of  broken  line)  obscure  the  comparison  not  only  with  the  other  computational 
data  (turbulence  models)  but  also  with  the  pressent  study.  Further  investigation  is 
required,  however,  to  determine  the  cause  of  these  deviations. 

Vorticity  variations  at  different  planes  are  shown  in  Figs.  9  through  11.  The 
contours  of  vorticity  component  in  the  streamwise  planes  (y-z  planes)  in  the  x-direction 
with  each  plane  identified  as  a,b,c,d,e  are  shown  in  Fig.9.  The  corresponding  velocity 
vectors  are  plotted  on  the  right-hand-side.  Clearly,  the  vortex  stretching  occurs  toward 
downstream  with  the  evidence  of  separation  shocks,  slip  lines,  and  vortex  centers  close  to 
the  wall.  These  physical  phenomena  become  more  significant  toward  downstream  in 
agreement  with  the  skematics  shown  in  Figs.fib  and  c. 

Fig.  10  shows  the  contours  of  vorticity  component  in  the  spanwise  vertical  planes 
(x-z  planes)  in  the  y-direction,  with  each  plane  idendified  as  a,b,c,d.  The  vortex  stretching 
occurs  again  toward  downstream  and  moving  upward  along  the  shock.  The  growth  of 
vorticity  is  concentrated  within  the  boundary  layer  close  to  the  wall. 

In  Fig.  1 1  the  spanwise  horizontal  plane  vorticity  contours  are  presented  at  various 
locations  (a  :  26g,  b  :  c  ;  0.5^;,)  where  S„  is  the  boundary  layer  thickness.  It  is  seen 

that  vorticity  increases  toward  the  wall  with  its  intensity  increasing  toward  downstream  as 
expected. 

m  CONCLUSIONS 

Based  on  the  preliminary  results  obtruned  for  the  direct  numerical  simulation  using 
the  MEI  Taylor-Galerkin  Legendre  polynomial  spectral  element  method,  it  appears  that 
our  original  goal  for  the  initial  attempt  has  been  successfully  achieved.  The  main  concern 
was  to  make  certain  that  the  irregular  node  cormectivity  of  the  h-p  process  can  numerically 
be  implemented.  Elaborate  data  structure  schemes  which  have  been  developed  are  the 
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major  factor  for  these  achievements.  In  addition,  the  preliminary  results  for  three- 
dimensional  computations  of  sharp  fin  shock  wave  turbulent  boundary  layer  interactions 
are  satisfactory. 

There  are  still  many  more  tasks  remaining  unexplored.  They  include:  (1) 
verification  of  3-D  fluctuations,  unsteadiness,  and  turbulent  micro-scales  as  related  to 
turbulent  Mach  number,  turbulent  Prandtl  number,  and  turbulent  Reynolds  number,  (2) 
characterization  of  compressibility  effects  and  relaminarization,  (3)  energy  spectrum  data 
versus  fi'equency  domain  and  complete  3-D  turbulent  statistics,  (4)  laminar-turbulence 
transition  instability,  (5)  reliable  optimal  control  of  h-p  interactions  with  Legendre 
polynomials,  (6)  temporal  and  spatial  dependency  of  implicitness  parameters,  among 
others.  They  constitute  challenging  future  tasks  in  years  to  come.  In  summary,  it  is 
concluded  that  the  direct  numerical  simulation  for  turbulent  compressible  flows  with  the 
Legendre  polynomial  spectral  element  method  appears  to  be  promising. 
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Appendix  A 

AMALOGY  BETWEEN  MEI  SCHEME  AND  FDM-T7D 


For  simplicity,  consider  a  one-dimensional  Euler  equation 

^  =  0  (Al, 

The  MEI  equation  for  (Al)  with  derivatives  written  in  terms  of 
nodes  u^,  Ui.^,  and  U1.2  becomes 


Auf*^ 

At 


Ax^ 


-  Auf.*!^)  > 


2AX^  ^ 


-  2AUiri 


(A2) 


where  C  is  the  Courant  number,  C  =  aAt/Ax. 

The  FDM-TVD  equation  for  (Al)  may  be  written  as 


du^ 

dt 


111 

Axi 


1 

2 


(Ui-  Ui.i) 


tUi.i- 


uj 


(u 


l‘l 


Ui) 


(A3) 


with 

a'  =  max  (0,a)  =  1/2  (a  +  |a|) 
a'  =  min  (0,a)  =  1/2  (a  -  jaj) 

Introducing  an  implicitness  parameter  s  for  the  time  derivative 
on  RHS  of  (A3)  in  the  form 


Ui  =  uf  +  sAu"*^ 


(A4) 
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Substituting  (A4)  into  (A3)  and  assuming  that 

a-  =  0,a*  =  a,Y*  i  *  Y*  3  =  Y 

i-1  1-1 

for  simplicity,  wa  obtain 

=  --S^(Aur"-  -  lIS4x(Aur"-  2Auf:^  Auf.\^) 

At  2Ax^  ^  ^  2Ax2  '  1  1  “X  1; 

-  fill)  -  -  2fi“i  ►  fi“2)  (A5) 

Ax  2Ax^ 

Comparing  (A2)  and  (A4)  reveals  that 

Si  *  -I'  S2  =  -  Y.  C  =  -7,  Sj  =  s  =  2S2. 

It  is  clear  that  the  MEI  formulation  and  FDM-TVD  scheme  are 
emalogous;  in  fact,  they  are  identical  under  the  assumptions  made 
above.  The  implicitness  parameters  s^  and  S2  in  the  MEI  scheme 
play  the  role  of  TVD  limiters,  7.  However,  the  implicitness 
parameters  s,  and  s,,  beyond  the  concept  of  TVD  scheme,  are 
expected  to  govern  complex  physical  phenomena  such  as  turbulent 
boundary  layer  interactions  with  shock  waves,  finite  rate 
chemistry,  widely  disparate  length  and  time  scales, 
compressibility  effects  in  high  Mach  number  flows,  etc. 
Undoubtedly,  all  implicitness  parameters  are  flowfield  dependent 
with  s.,  and  s,  associated  with  Mach  numbers,  s,  and  s^  with 
Reynolds  numbers  and  Damkohler  numbers,  between  upstream  and 
downstream  nodes  within  finite  elements. 
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Fig.  3  Perturbation  velocity  components  for  the  flat 
plate,  page  37,  1993  progress  report. 
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Fig.  4  Reynolds  stree  components  for  the  flat  plate 
page  37,  1993  progress  report. 
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Pig.  5  Turbulent  kinetic  energy  for  the  flat  plate 
page  37,  1993  progress  report. 
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streamline  vorticity  contours  and  velocity  vectors.  ( t=0.3965ms) 
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